Surface Tension in Unitary Fermi Gases with Population Imbalance 
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We study the effects of surface tension between normal and superfluid regions of a trapped Fermi 
gas at unitarity. We find that surface tension causes notable distortions in the shape of large aspect 
ratio clouds. Including these distortions in our theories resolves many of the apparent discrepancies 
among different experiments and between theory and experiments. 



Experimentalists are now using dilute gases to control- 
lably study the properties of strongly interacting systems 
of superfluid fermionic atoms |lj, [2J, y, |4|. Recent ex- 
periments have examined the exotic circumstance where 
atoms with two different hyperfine spins [denoted up and 
down] are placed in a harmonic trap, but the number of 
spin-up atoms, N_t is greater than the number of down- 
spin atoms 7Vj [2l- Li Li Spin relaxation is negligible 
in these experiments, so over the entire time of the ex- 
periment, the system is constrained to have a fixed po- 
larization P=(7V T - JV;)/(7Vf + JV;). Understanding the 
structure of (s-wave) superfluidity in this polarized en- 
vironment is an important endeavor with a long history 
[a, |6|, III l| and direct relevance to neutron stars, thin-film 
superconductors, and color superconductivity. In this 
paper we use the concept of surface tension to quantita- 
tively explain controversial features seen in the density 
profiles of strongly interacting trapped polarized Fermi 
gases HHQ. 

The simplest theories of trapped Fermi gases [{J, Hfl 
llll Il2l Il3| (most relying on local density approximations 
[LDA] and assuming zero temperature) predict that the 
atomic cloud phase separates into a central superfluid re- 
gion, in which the density of both spin species are equal, 
surrounded by a polarized normal shell 14]. This ba- 
sic structure was observed in two separate experiments 
[2J, Li L! , however some experimental details are at odds 
with these theoretical predictions. For P > 0.1, the Rice 
experiments |4| find a double peaked axial density dif- 
ference, nf'(z) = J dxdy [n-f(r) — nj.(r)], where n-j-n(r) 
is the density of up and down spin atoms. In a previ- 
ous paper |2(, we argued that this structure pointed to 
a breakdown of the local density approximation, despite 
the fact that dimensional arguments suggested that the 
LDA should work well. Conversely, the results of the 
MIT experiments |2], y] are fully consistent with a local 
density approximation, but show a polarization driven 
superfluid-normal phase transition at P ~ 0.70. This 
phase transition was not seen in the Rice experiments and 
is not found in most theories at unitarity jjj, |l(J, llll Il2j ] . 
Here we show that surface tension in the boundary be- 
tween normal and superfluid regions distorts the cloud 
in exactly the right way to account for the unusual fea- 
tures seen at Rice. We also show that surface tension 
plays a much smaller role in the MIT experiments, where 



the atomic clouds are larger and more spherical, and we 
are thus able to account for the fact that the MIT ex- 
periment is consistent with the local density approxima- 
tion. Finally, we show that for P > 0.7, the Rice data 
shows a sudden drop in surface tension. Since such a drop 
would be expected if the system underwent a superfluid- 
normal phase transition, this observation may reconcile 
the apparent differences in the experiments. We cur- 
rently lack a quantitative theory of the superfluid-normal 
phase transition at unitarity. 

In this letter we consider the unitary regime, where 
the scattering length is infinite and the only lengthscale 
in the problem is the interparticle spacing. Taking a two- 
shell structure, with a superfluid core and a normal fluid 
shell, we model the free energy of a trapped gas as 
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where J s / N represents the integral over the super- 
fluid/normal region, f g corresponds to an integral over 
the boundary, f g i„ = — Jn^'^'dfi represent the free 
energy density of the superfluid/normal gas and a rep- 
resents the surface tension in the boundary. The energy 
densities are a function of the local chemical potentials 
ix[f) = \[i X ir) +/i;(r)]/2 = /i - V(r), and h = bj T (r) 



Ml(?~)]/2, where V(r) 
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is the trapping potential, with A sa 50 for the Rice exper- 
iments and A w 5 at MIT. The shape of the boundary, 
and the parameters /io and h, are determined by min- 
imizing eq. (Q with respect to the boundary with the 
constraint that Ny^ — J s d 3 7 



,(») 
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approach is a generalization of one used by Chevy [lfj, 
where the boundary term was absent. Universality allows 
us to write the free energy density as 
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where ( s = 1/(1 + Pf' 2 , ( n = 1/2, and j3 w -0.545 
is a universal many-body parameter [l5j. The rele- 
vant chemical potentials are p s (r) — /j,(r) and /U„(r) = 
fi^(r) = fj,(r) + h. The density of each spin component 
is w\,i = —df/dn. The fact that the particle spacing 
is the only length scale constrains the surface tension to 




FIG. 1: TOP: Comparison of the minority and majority 
components radii with the Rice experiments 4]. Squares 
and crosses are the majority and minority components ob- 
served in U, determined by looking for where the density 
of each component vanishes. The solid line is our theoret- 
ical prediction of the radii with zero surface tension, while 
the dashed line includes finite surface tension. The radii 
are normalized with a non-interacting Thomas-Fermi radius, 
Rtf — \/^f/bz), where e/ = TuJJftN) 1 ' 3 with average trap 
frequency Zu = {tu\uj z ) iri and N - (iV T + N{]/2. Note, both 
the normalization and the fitting procedure differ from the 
one used in figure 3 of Q. 



have the form a = (h /2m)n s (r)g(5P/ 'P), where g is 
a function of 6P, the pressure discontinuity across the 
domain wall, and P, the pressure on the superfluid side 
of the domain wall and n s is the density on the super- 
fluid side. Introducing a universal numerical parameter 
?7, we approximate g by its value at zero pressure drop, 
q(0 ) = 77, based on estimating that SP/P < 1.8 x 1CT 3 

We determine r\ in two ways. First, as detailed in the 
appendix, a mean-field theory gradient expansion yields 
j] w 0.9 x 10~ 3 . Second, we use a fitting scheme where 
we minimize Eq. ( ^) for a series of candidate 77's. We 
find that r\ — 1.0 x 10 -3 matches the Rice group's exper- 
imental data for the axial density difference at P=0.53. 
Given the uncontrolled nature of the mean field approx- 
imation, we believe that the similarity of the two results 
is purely coincidental. We use r\ = 1.0 x 10~ 3 for all of 
our predictions. 

To simplify the minimization of Eq. (Q with respect 
to the boundary, we make the ansatz that the boundary 
is an ellipsoid with semi major and minor axes ~p and z. 
Within this ansatz we analytically calculate the free en- 
ergy [for brevity we omit the expressions] . We minimize 
this expression with respect to the parameters ~p and z. 

To estimate the distortions, one expands Eq. Q for 
small distortions: ~p = p (l + S p ), and ~z — Zo(l + 5 Z ), 
where po and zq are the lengths of the axes in the absence 
of surface tension. We take S p and 5 Z to be order of S. Di- 
mensional analysis gives, SF/(h 2 /2m) ~ A rjn 4 ^ 3 zopoS + 
B n 5 / 3 zopQ5 2 , where A and B are constants. Assuming 
that po scales with the radial Thomas-Fermi radius, the 
size of the distortion is then 5 ~ l/(p n 1 / 3 ) ~ (A/iV) 1 / 3 . 

FIG. n shows the calculated axial radii as a function 
of polarization. We compare our predictions to radii that 




FIG. 2: LEFT: Axial density difference n ( f' \z) = 

2n J dpp\n-j(z,p) — ni(z,p)] of zero temperature harmoni- 
cally trapped unitary Fermi gas in units of [10 6 cm ]. Fig- 
ures (a), (b), and (c) represent polarization P = 0.14, 0.53, 
and 0.72 respectively. The grey points are the experimental 
data from reference H. The P = 0.14 and P — 0.53 data 
previously appeared in figure 2 of |4J. The P — 0.72 data 
corresponds to one of the points in figure 3 of ref 4J. The 
dashed line is the zero surface tension density, while solid line 
is the finite surface tension density using 77 = 1.0 x 10 -3 . 
RIGHT: Comparison of the axial density of minority compo- 
nent n^ (z) — 2ty j dp pni(z, p). 



we extract from the data used in FIG. 3 of ref. [4|. We 
extract the radii by fitting the wings of the axial den- 
sity distributions to a piecewise linear function of the 
form n(z) = w (1 — z/R)8(R — z), where 6(x) = 1 for 
x > and otherwise, and w and R are fitting parame- 
ters. This fitting procedure lets us accurately determine 
the edge of the cloud, while the radii extracted in [4| 
correspond to an average radius, and are systematically 
larger those extracted by our method. We see that for 
0.1 < P ?S 0.7 the experimental data is in excellent agree- 
ment with the finite surface tension theory. For P > 0.7, 
the data used for FIG.n a PP ears to be inconsistent with 
rj = 1.0 x 10~ 3 . We speculate that the deviation may be 
due to the superfluid-normal phase transition observed 
at MIT |2j, |3j . We caution, however ,that there are large 
fluctuations in the radii seen in ref. |4| especially at large 
P. More work is needed before definitive statements can 
be made. The disagreement of the radii below P c « 0.1 
is probably attributable to finite temperature effects [lj| . 
In FIG. [3and|21wc compare our predicted axial den- 
sity profiles with representative data from ref. |^|. As 
demonstrate in the left panel of FIG.[3 the finite surface 




FIG. 3: Comparison of the axial majority component density 
ni a (z) = 2ty J dp pn^(z, p) in units of [10 6 cm -1 ] with exper- 
imental data of reference Q. Figure (a) and (b) represent 
polarization P — 0.53 and P — 0.72 respectively. Symbols 
carry the same meanings as in Fig [5] Notice that the solid 
line and the dashed line coincide, indicating surface tension 
has no effect on the majority densities. 



domain wall energy in this approximation is given by, 
E ( l ) = S% dx b I d * A I 2 +E bcs (A,h,n)-E n (h,n)], 
where I, the size of the domain wall, will be determined 
variationally, along with A(x), the superfluid order pa- 
rameter. The energies in the superfluid and normal 
phases are given by jjj, 



E bcs (A, /j, M ) = — 2 / k 2 dk{~h + E k ) (3) 
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tension theory captures the observed double peak struc- 
ture in axial density difference for P < 0.7. The only free 
parameter in this calculation is 77, which as previously de- 
scribed we set by fitting to the P = 0.53 data. A close 
examination of FIG. \^c) reveals that the P=0.72 data is 
not fit quantitatively by either the finite surface tension 
or zero surface tension theory. As previously discussed, 
we suspect that the central region may not be superfluid. 

As illustrated in FIG. |3 surface tension has almost 
no effect on the axial density of the majority component. 
The smallness of the effect is to be expected because the 
discontinuity in nj at the domain wall is much smaller 
than the discontinuity in rij.. Alternate explanations of 
the double-peaked axial density difference, such as an- 
harmonicities [lj|, would cause distortions in n-f instead 
of n j, and are not completely consistent with the experi- 
mental data |19j . 

We also calculated, but do not show here, density pro- 
files for parameters corresponding to the MIT experi- 
ments. We find that surface tension has a negligible ef- 
fect on the density profile, consistent with the fact that 
(A/TV) 1 / 3 is 10 times smaller than at Rice. 

We wish to emphasize how surprises it is to see surface 
tension, a phenomenon generally associated with liquid 
in a gas. This observation opens the possibility of other 
surface tension related effects in cold atoms. In particular 
the surface tension could have a large effect on collective 
modes and expansion. We speculate that surface tension 
should play a role in the physics of analogous systems, 
such as nuclear matter at high densities and quark-gluon 
plasmas. 

This work was supported by NSF grant PHY-0456261, 
and by the Alfred P. Sloan Foundation. We are grateful 
to R. Hulet and W. Li for very enlightening discussions, 
critical comments, and for providing us with their exper- 
imental data. We also thank M. Zwierlein for providing 
us with their latest results [3j, and for insightful critical 
comments. We thank T.-L. Ho for critical comments. 

APPENDIX: In this appendix we calculate the do- 
main wall energy at the superfluid-normal interface by 
applying a gradient expansion to mean field theory. The 



E n (h,n) 
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2m 
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where k±(r) = (±\/h 2 - A 2 + /1) 1 / 2 , E 2 = e 2 + A 2 , e fc = 
h k 2 /2m — fj, and, a s is the s-wave scattering length. 

In order to calculate the coefficient of the gra- 
dient term 7, we begin with the action S — 
So + Sinti where the free fermions action is So = 
Y,<, Jo drf d 3 FxPUd T - av + h 2 V 2 /2m]A*, the in- 

teraction is Si n t = — U f dr f d 3 fif)l.^il>hip r iip r -[, the 
atomic Fermi fields are i/w, imaginary time is r, 
the inverse temperature is (3 = \/T and the at- 
tractive interaction between Fermi atoms is —U with 
U > 0. After the usual Hubbard-Stratonovich de- 
coupling of the interaction term [20j and integrat- 
ing out the Fermi fields, the partition function is 
written as Z = J DADA*exp[-S ef f(A, A*)], where 
S eff (A 1 A*) = j: q JA( q ,u n )\A( q )\ 2 + B( q ,u J )\A( q )\* + 

...] and A( 9)Wn ) = (1/U -TJ2 n fd 3 k/(2n) 3 G ] -(k + 
q /2,L0 m )Gi(-k + q/2,u> m + w n )) with G a (k,ujy l = 
iuj — h k 2 /2m + \x„. We assume that the dominant mo- 
mentum dependence comes from the term which is lowest 
order in superfluid order parameter. We sum over Mat- 
subara frequencies, defining A{ q ) = *52 n A( q ,u} n ) with 
w n = (2n+ 1)ttT. In order to suppress the ultraviolet di- 
vergences in the theory, we regularize |2l| the interaction 
with the s-wave scattering length by 1/U = m/ATrh 2 a s + 
d 3 k/(2Tr) 3 m/h k 2 . We then expand A(q) to second order 
in q and take the zero temperature limit, finding A( q ) = 
m/(4irh 2 a s ) - m^/JI/(4TTh 2 ) + m q 2 / \Z2ttTi 2 Jji) + C(g 4 ), 
which means that 7 = m/(32irh y/jT). 

Taking the ansatz A(x) = (Ao/2)[tanh(x/^) + l], where 
Ao is the value of A on superfluid side of the domain wall, 
we numerically minimize the surface energy E(l) with re- 
spect to / to find the size of the domain wall l m and the 
domain wall energy. We find l m < kj , supporting our 
treatment of the domain wall as very thin, but calling 
into question the validity of our gradient expansion. Note 
that we do not expand Eb cs in powers of A, but work with 
the exact expansion. Since we are considering a domain 
wall between a region where A = and A w O(Ef), 
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FIG. 4: The value of 7/ as a function of (fc/a s ) _1 in zero 
temperature BCS approximation. At unitarity r\ = 0.0009 
is independent of the density. The Fermi wave vector kr is 



any expansion in A would require going to high order. 
To even capture the topology of the free energy surface, 
one must expand to sixth order. Thus, prev ious calcula- 
tion of surface tension, such as Caldas's |22( recent work, 
which are based on fourth order expansions are not rel- 
evant to the physics described here. By repeating this 
calculation at different a s , and solving the BCS num- 
ber equation and gap equation |9J, we find the quantity 
7] = 2mE(l m )/(fi n s ) as a function of a s and n s , where 
n s is the density on the superfluid-normal interface. In 
the limit of a s — > 00, we find 77 = 0.9 x 10~ 3 , indepen- 
dent of the density and polarization. However, as seen in 
FIG. QJ 77 has density dependence away from unitarity. 
As a s — > + , 77 grows larger, hence the effects of surface 
tension is stronger. Therefore, in the strong BEC regime, 
domain walls become energetically prohibitive and the 
phase separated atomic system is unstable against phase 
coexistence Lillsl lllj . Recent theoretical work by Imam- 
bekov et al |2j| studied the role of gradient terms in this 
deep BEC limit. 

The value of r\ obtained from fitting to experimental 
data agrees well with our mean field calculation. We 
believe that this agreement is coincidental as mean field 
theory is not expected to work well at unitarity. We 
also note that the experiment is performed slightly away 
from the resonance where the mean field approximation 
predicts a weak density dependance of r\ [24j . 
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